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Abstract 

Recently, Papadakis et al. na proposed an efficient primal-dual algorithm 
for solving the dynamic optimal transport problem with quadratic ground 
cost and measures having densities with respect to the Lebesgue measure. 
It is based on the fluid mechanics formulation by Benamou and Brenier [2 
and proximal splitting schemes. In this paper we extend the framework 
to color image processing. We show how the transportation problem for 
RGB color images can be tackled by prescribing periodic boundary con¬ 
ditions in the color dimension. This requires the solution of a 4D Poisson 
equation with mixed Neumann and periodic boundary conditions in each 
iteration step of the algorithm. This 4D Poisson equation can be effi¬ 
ciently handled by fast Fourier and Cosine transforms. Furthermore, we 
sketch how the same idea can be used in a modified way to transport pe¬ 
riodic ID data such as the histogram of cyclic hue components of images. 
We discuss the existence and uniqueness of a minimizer of the associated 
energy functional. Numerical examples illustrate the meaningfulness of 
our approach. 


1 Introduction 

Recently, methods from optimal transport (OT) have gained a lot of interest 
in image processing. In one of the first applications, the Wasserstein distance 
(earth mover distance) has been successfully used for image retrieval [19 and 
since then it has been applied to many other tasks as color transfer mm, 
(co)segmentation [TO} [T4l [20] . the synthesis and mixing of stationary Gaussian 
textures pQ and the computation of barycenters Hsus]. 

The basic problem, going back to Monge (1746-1818), can be formulated as 
follows: Given two probability spaces (T,/io) and (3^,/ii) and a nonnegative 
cost function c(x, y) on A x y, find a transport map T: X —>• y that transports 
the mass of fi o to the mass of yi at minimal cost, i.e., T minimizes 

/ c(x, T(x)) dyo(x) subject to /io°T _1 =/ii. (1) 

J x 

One of the major limitations in applications using OT is the fact that it is in 
general not known whether a solution of problem 0 exists and even in this 
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case the computation of the optimal map T is usually a demanding task (except 
very few cases, e.g. OT on R with convex costs). We focus in the following on a 
specific instance of problem ([!]), namely if X — y = R d , c(x, y ) = \ \x — y\ 2 and 
/io and fi i are absolutely continuous w.r.t. the Lebesgue measure, that means 
there exist probability density functions f° and / 1 with 

IH(A)= [ f*{x) dx, z = 0,1, AeB(R d ). 

J A 

In this case there exists a unique optimal transport map T that transports f° 
to Z 1 , see, e.g., li- Instead of considering a time independent, “static” 
mass transportation problem one may alternatively consider the geodesic path 
between the two measures w.r.t. to the Wasserstein metric (the so called dis¬ 
placement interpolation w- While the static problem can be seen as a distance 
problem (find the minimal distance between the probability measures /io and 
/ii), the dynamic problem can be interpreted as a geodesic problem (find an 
optimal path between /io and /ii). For the L 2 ground cost this geodesic is ob¬ 
tained by linear interpolation between the identity and the optimal transport 
map T, i.e., /it = /io ° T t -1 , where T t = (1 — t) Id +tT. Benamou and Brenier j2] 
gave the following equivalent formulation of the dynamic OT problem in terms 
of fluid mechanics: minimize 

[ [ bOM) \v(x,t)\ 2 dx dt, (2) 

subject to Ut€[o,i] SU PP/(V) bounded, /(•, 0) = /°,/(-, 1) = / 1 and d t f + 
diy x ifv) = 0, (f,v) sufficiently smooth. Substituting m = fv, this problem 
becomes convex and can be treated by respective algorithms. 

In addition to the Dirichlet boundary condition for the time interval, problem © 
needs to be equipped with spatial boundary conditions in practical applications 
(appearing in the momentum variable m) . A natural choice which was also used 
in m for gray-value images are Neumann boundary conditions. In this paper 
we want to deal with color RGB (red, green, blue) images. More precisely, we 
consider a m x n RGB image as a 3D object of size mxnx 3 with the color values 
in the third dimension, i.e., we interpret these images as (realization of) a 3D 
density function. To use Neumann boundary conditions for the color dimension 
is certainly not a good idea since the solution can depend on the ordering of the 
color channels. Therefore, we suggest to establish periodic boundary conditions 
in the third dimension. Note that we have Dirichlet boundary conditions in 
time, and Neumann plus periodic spatial boundary conditions. Fig. [l] illustrates 
the effect of the different boundary conditions. In Section [2] we explain the 
discretization of problem Q emphasizing the mixed boundary conditions. In 
particular we deal with the existence and uniqueness of minimizers. Similarly 
as suggested in [12 we apply a primal dual algorithm to find a minimizer in 
Section [3j In contrast to pjj each iteration step of this algorithm requires the 
solution of a 4D Poisson problem with mixed Neumann and periodic boundary 
conditions which can be efficiently computed by applying FFTs and fast cosine 
transforms. In Section [4] we apply our findings to the dynamic OT of RGB 
images which was the initial motivation of this work. Another application is 
given in Section |5j Here, the (cyclic) OT is applied to HSV (hue, saturation, 
value) images, where only the cyclic hue component is transported. For further 
details we refer to 0. 
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Figure 1: OT results for transporting a red Gaussian into a blue one with different 
boundary conditions for the third (color) dimension. First row: Neumann 
boundary conditions; the color is transported from red over green to blue. 
This changes if we permute the RGB channels. Second row: periodic 
boundary conditions; the color is transported from red over violet to blue 
which is more intuitive and does not change for permuted color channels. 


2 Model for dynamic OT with mixed bound¬ 
aries 

Rewriting © for m = fv, we see that the geodesic path between measures with 
probability densities f° = /(•, 0) and / 1 = /(•, 1) has density /(•,£) fulfilling 

argmin / / J(m,/)dxd£, (3) 

(mj)ecJ[ 0,1] J[ 0,l] d 


where 




0 if = (0j,0), 

Too otherwise, 


c ■■= {( f,rn ) : d t f + div x m = 0,/(-,0) = f°, /(•, 1) = f 1 } 

with appropriate boundary conditions. In the following, we describe the dis¬ 
cretization of problem © for one spatial dimension with cyclic spatial boundary 
conditions. The discretization of the continuity equation demands the evalua¬ 
tion of discrete partial derivatives in time as well as in space. In order to 
avoid solutions suffering from the well-known checkerboard-effect (see for in¬ 
stance [T^) we adopt the idea of a staggered grid as in m ■ 

Discretization (Spatial ID): We consider the values of / at spatial cell mid¬ 
points , j = 1,..., n and time k = 1,... ,p — 1, and denote the corre¬ 

sponding array by (/(j — |, k G M n,p_1 . The boundary values are as¬ 
sumed to be fixed in f° := (^f C~^ 2 ,0)^ and / 1 := (^f C~^ 2 ,1)^ , where 

P > 0, i = 0,1 and ||/°||i = ||/ x ||i- We can skip the normalization ||/°||i = 1 
here. The values of m are taken at the cell faces j = k, ..., n — 1 and time 

k ~p / 2 , k = 1,... ,p and we consider the array (ra(j, fc — |)) *’£ =1 G 

where ft = 1 for Neumann boundary conditions and ft = 0 for periodic boundary 
conditions, see Fig. |2j To give a sound matrix-vector notation of the discrete 
minimization problem we reorder m and / columnwise into vectors / = vec(/) G 
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M n ( p —i), m = vec(m) E R( n ~ K ) p . Since it becomes clear from the context if we 
deal with arrays or vectors we use the same notation. The derivatives in C are 
approximated by forward differences and the integral in ©> by a midpoint rule, 
where the midpoints (u(j — \,k— |))J’^ =1 and (v(j — |, k — are com¬ 

puted by averaging the neighboring two values. This results in the following 



Figure 2: Discretization grid for the dynamic OT problem for ID signals /: horizon¬ 
tal direction for time; vertical direction for space; • given boundary sam¬ 
pling nodes f° and f 1 , ■ sampling nodes for f((j — 1/2)/n, k), j m 1,..., n, 
k — 1 ,.. . ,p — 1 , □ sampling nodes for m, ’+’ quadrature nodes. 

discrete model: 


argmin \\J(u,v) ||i, 

(4) 

m,f,u,v 


subject to Sm'>tT' = u, Sf/ + =v i 


(Av4Df)(7) =/4 

(5) 


A 


We denote 

Cd ^ { (?) G ® (n “ K)P+n(j,_1) = (Dm\D f ) (j) m 
The involved vectors are defined as 

fb ■= ^((/ 0 ) t , o „ (p _ 2 ),(/ 1 ) t ) t , 
fb ~p{{f°TAn{p-2),-U X Y) T 

and the matrices using the Kronecker product 0 as 


S F ~S*®I n , D F ~-Dl®I n , 


Sm •= 


I p 0 Neumann, 

4 ® Sl tPer periodic, 


and 


Sp ■— 2 


Dm ■= 


f\ ; 


f 1 \ i 


/ 4' 
l 4' 


• —DT Neumann, 


D 


n,per 


periodic, 


1 \ 


110 
i i ) 


e 


i i / 
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D n 


0 
-1 


D p := p 


V 


-1 10 
-1 1 / 


Since we have to be slightly careful concerning the uniqueness of the solution 
in the periodic setting we provide the following proposition. 

Theorem 1. The discrete dynamic transport model 0 has a solution. 

Proof. For periodic boundary conditions and even n, we have AT(S'm) = {w<g)l n : 
w G with l n := (1, —1,..., 1, —1) T G M n , and A/"(-Sm) = {0 np } otherwise. 
The constraints in Cd can be rewritten as 

D F f = fb ~ D Mm , / = D ] F f~ - D ] F D M m 

with the Moore-Penrose inverse D F = (D] p Df)~ 1 D f . Then we obtain 

argmin || J(u, = argmin || J (X(ra), Y (m))^ 


withX(ra) := Y(m) := — SFDpDMm+SFDpf^+f^ . Let ||m||2 —^ +oc 

and assume that \\J (X(m),Y(m)) ||i is bounded. Then each of the quotients 
is bounded, i.e., there exists c > 0 such that X(m)f < c\Y(m)i\ and therefore 
||X(m )||2 < c\\Y(m)\\i. Thus, in the case J\f(S m ) = {0 np }, we get 

(VII^||I)|H||<||X(m)|||<c||F(m)|| 1 (6) 

< c\\S f D^ f Dm IIi ll^ll i +c, 

which is not possible as ||m ||2 Too. In the case A/"(5 m) = {w ® l n • w E 
M p } we use the orthogonal splitting m = uir + w <g> l n , where tur G 7 Z{S r f /[ ). 
Straightforward computation shows that SeD^Dm^^^-u) = w® l n with some 
w G M p . Since Y(m) > 0, we conclude 

|| - SfD^Dm^R + + / 6 + ||oo > \\™ ® f n 11 oo • (7) 

If mR remains finite, then In ® w remains finite as well. Since the kernel of 
SfDfD m consists of constant vectors, this is only possible if re is a multiple of 
1 p. But in this case J p (X(m),Y(m)) has a finite value which is reached. It 
remains to consider the case ||ra #||2 —>• +oc. Then we obtain similarly as in ^ 
that 


(VII^Lk(s-)III)l|mfillI < II^MIII < cIVMIIi 

< c|| - S F D ] F D M m R + SpD^F + fb Hi + c \\™ ® || 1 - 

By 0 we see that this is not possible as ||m ^||2 —^ +oc. In summary, we 
have that ||J (X(m), Y(m)) ||i is coercive and since it is also proper and lower 
semi-continuous, it has a minimizer. □ 
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Unfortunately, J(ra, /) is not strongly convex on its domain. As it can be seen 
in the following lemma it is even not strictly convex. 


Lemma 1. For any two minimizers (mi,fi), i — 1,2 the relation 

Sm^i Sm^ 2 


S F fi + f b S F f 2 + f b 


holds true. 

Proof. The function J{u,v ) is the perspective function of the strictly convex 
function 'ip = | • | 2 , i.e., J{u , v ) = vip (^), see, e.g., 0. For A G (0,1) and ( Ui,Vi ) 
with ^ > 0, i m 1, 2, we have (componentwise) 


Xu\ + (1 — X)u 2 


J (X(ui,vi) + (1 - X)(u 2 ,v 2 )) = (Xv! + (1 - X)v 2 ) ip 


Xv 1 + (1 - X)v 2 



and if — ^ — by the strict convexity of ip, 

V\ ' V2 J J “ ‘ 

J (A(«i, ui) + (1 - A)(it 2 ,i; 2 )) < AJ(«i,vi) + (1 - X)J(u 2 ,v 2 ), 


which proves the assertion. 


□ 


Remark 1 . For periodic boundary conditions, even n and f 1 = f° + 7 I n , 
7 G [0,min/°) t/ze minimizer of Q zs not unique which can be seen as follows: 
obviously, we would have a minimizer (m, /) if m = w (g) l n £ -XP(Sm) f or 
some w G and there exits f > 0 which fulfills the constraints. Setting 
jk/p ._ _ i/2,k)f =1 , k = 0, -..,p, these constraints read —2 pw <S> l n = 

p(f(k~l)/P _ fk/p)P =im Thus; 

any to G l p stzc/z that 


f 1/p = f° + 2wil„, / 2/p = f° + 2(u>i + W 2 )l n , • • • , 
f 1 = /° + 2(wi + w 2 + ... + 


are nonnegative vectors provides a minimizer of ©• ITe conjecture that the 
solution is unique in all other cases, but have no proof so far. 


3 Minimization Algorithm 

We apply the primal-dual algorithm of Chambolle and Pock [4] in the form of 
Algorithm 8 in [3]. Step 1 requires the projection of 



onto Cd which is given by 


n c d (a) = a- A T (AA T )\Aa - f b ), 
(AA T ) t = Qdiag(A : ,-)(2 T , 
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Algorithm 1: PDHG Algorithm for solving Q 

Initialization: ra (0) = 0 np , / (0) = 0 n(p _i), b$ = bf^ = b$ = 5^ 0) = 0 np; 
0 = 1, T, <T with T<7 < 1. 

Iteration: For r = 0,1,... iterate 


1 . 


W r+1 >\ 

j(r+!) ) 


2 . 


0+l)\ 

0+1) J 


3. & +1) 
fc (r+1) 

4. b£ +1) 
fc (r+1) 


:= argrnm — 

( m,f)ec d 




+ TCT 



:= argmin||J(tz,v)||x +^|| 

u,v 



/5m 0 \ /m( r+1 )\ _ ( 0 

VO 5f/ V/( r+1) J [f~ 

:= bf^ + SM , m ( ' r+1 ' 1 — w (r+1) 
:= 6 V r) + ^/ (r+1) + /fc + -« 
:= 6 ^ + 1 ) +^ + 1 ) - 6 ^) 
:= # +1) + 0 (b}T +1) - b^) 



O+i) 





where AA T has the spectral decomposition AA T = Q diag(Aj)Q T and A j := 1/A j 
if A j > 0 and zero otherwise. For A in <§ the application of (AA T ) 1[ requires 
the solution of a 2D Poisson equation. In case of periodic boundary conditions 
we get 


where 


Since 


AA T = I p ® D£ jper D niPer + DpDp ® I n 
— Ip ® n A n? p er H - p /Ap ® / n , 


( 2 - 1 

-1 2 -1 


V -1 

/ 1 - 1 
-1 2-1 




1 - 


V 0 


-i 


-1 2 -1 
-1 2 ) 


-1 2 -1 
-1 1 / 


A n,per 


— F n diag(g n?per )F n , Qn,per •— I 4 sin 


• 2 


n 


n — 1 


3 = 0 


with the Fourier matrix F n := (e 27rljfc )^ fe ^ 0 an( 4 


A p = Cp diag (q p )C p , q p := I 4 sin 2 ^ 


P -1 


with the DCT-II matrix 


Cr) := \ - ( Co cos 


j(2k + 1)7T 
" 2 p 


p-i 


j,/e =0 


-{1 


2 Pij=o 

1/V2 if j = 0, 

otherwise 
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t = 


we obtain 

AA T =(Cp ® t F n )(n 2 I p <g> diag(g n)Per ) 

+ p 2 (diag((/ p ) <E> I n )(C p ® F„), 

so that its pseudo-inverse can be computed by the FFT and the fast cosine 
transform. 


Remark 2. For the transport of general 2D RGB images we have analogously to 
solve a 4D Poisson equation with Neumann boundary conditions and a periodic 
boundary condition for the color channels. 


Step 2 of the algorithm can be computed componentwise as proposed in m- 
Setting a rn := SM^ r+V) + < 2 / •= 5V/( r+1 ) + b^p + f^~ we have to find 

componentwise 

2 

. U (7 . ® / \2 

argmrn - + -{u-a m ) + ^{v - a f ) . 

u,v z z 

Setting the gradient to zero yields 

u 1 v? 

- + a(u - a m ) = 0, - - -2 + &(v - a f ) = 0. 

V z v z 


Thus, 

crvam 
u = -- 

1 + (TV 

and v is the solution of the third order equation 


f(v) = 2(1 + av) 2 (v - af) - aa 2 m = 0 . 


This can be solved by few Newton steps which can be computed simultaneously 
for all components. Alternatively, one may use Cardan’s formula. 


8 









Figure 4: OT between hue histograms and histogram speciflcation at time t = j = 
0 . 8 . 

4 RGB Image Transport 

In our first experiment we consider the periodic color OT between two RGB 
images uo and u\ that are used as densities f° and / 1 respectively. At this 
point it is important to choose image pairs which have approximately the same 
mass (i.e. the overall sum of all pixels and color channels) as one needs to 
rescale the images such that 11 y* 0 11 x = H/ 1 ^. The results of four different 
examples are shown in Fig. [3j The first row shows an artificial example of 
the transport between one red Gaussian into a blue and a green Gaussian with 
smaller variance. In the second row, two polar lights of different color and shape 
are transported into each other. The third row illustrates how a topographic 
map of Europe is transported into a satellite image of Europe at night. Finally, 
the last row displays the transport between two cranes in Hamburg. All images 
are of size 100 x 100 x 3 and in each case, we used P = 32 time steps and 2000 
iterations in our algorithm. Note that, however, already after 200 iterations 
there are no visible changes any longer. The images are depicted at intermediate 
times t = |, i = 0,..., 8 . In all cases one nicely sees a continuous change of 
color and shape during the transport. 


5 Hue Histogram Transport 

In this section we perform color OT in the HSV space. We assume the final 
and target image have the same saturation and value such that only the cyclic 
hue component has to be transported: Assume we are given two images Ui , 
i = 0,1 which differ only in the hue component represented by their normalized 
histograms h l as empirical densities f\ i = 0,1. As the hue component is 
periodic, this fits into our setting. The intermediate histograms ft,*, t G (0,1) are 
then used to obtain the hue images via histogram specification. Together with 
the original saturation and intensity they yield the images ut , t £ (0,1). For the 
histogram specification of periodic data we have applied the analysis in m and 
the exact histogram specification method for real-valued data proposed, e.g., 
in n, Fig.ra shows an example, where the histogram of the hue component of 
a yellow flower is transported into the one of a red flower. The color changes 
gradually and in a realistic way which would not be the case if the periodicity 
of the hue histogram is not taken into account. 

Acknowledgement: Funding by the DFG within the Research Training Group 
1932 is gratefully acknowledged. 

*A11 images from Wikimedia Commons: AGOModra_aurora.jpg by Comenius University 
under CC BY SA 3.0, Aurora-borealis_andoya.jpg by M. Buschmann under CC BY 3.0, Eu¬ 
rope-satellite-orthographic, jpg and Earthlights_2002.jpg by NASA, Kohlbrandbriicke5478.jpg 
by G. Ries under CC BY SA 2.5, K6hlbrandbriicke.jpg by HafenCityl under CC BY 3.0. 
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